%%
data = load('C:\Manoj\projects\cej\stefan\EEJ_parameters_V3.0\EEJ_parameters_V3.0\CHAMP_scalar.txt');
lt1 = champ_lt_new(data(:,1),data(:,2));

ncount=1;
for i = 0.1:0.01:5,
    
    power_eej1(ncount,:) = inversion(data(:,1),data(:,6),i);
    nperiod1(ncount) = i;
        ncount=ncount+1;

end;

plot(nperiod1*24,sqrt(power_eej1(:,2).^2+power_eej1(:,3).^2));
hold on;

L = lt1 >9 & lt1 < 13 & data(:,10) > 3;
L = data(:,10) > 3;
clear nperiod2 power_eej2;
ncount=1;
for i = 0.1:0.01:5,
    
    power_eej2(ncount,:) = inversion(data(L,1),data(L,6)./1.9e-6,i);
    nperiod2(ncount) = i;
        ncount=ncount+1;

end;

plot(nperiod2*24,sqrt(power_eej2(:,2).^2+power_eej2(:,3).^2),'k');

% st=1;
% ncount=1;
% for i = 1:40000,
%     
%     if sum(L(st:st+127)) == 128,
%         
%         dataset(ncount) = st;
%         st = st+128;
%         ncount=ncount+1;
%     else,
%         st=st+1;
%     end;
% end;
%%


load c:\manoj\projects\ace\OMNI_ELEC_new.mat;
L = isnan(ace_all(:,2));
new_ace_data = interp1(ace_all(~L,1),ace_all(~L,2),ace_all(:,1));
new_ace_data_inter = interp1(ace_all(:,1),new_ace_data,data(:,1));


st=1;
for i =1:length(dataset),
    ace_seg(st:st+127) = new_ace_data_inter(dataset(i):dataset(i)+127);
    champ_seg(st:st+127) = data(dataset(i):dataset(i)+127,6);
    time_seg(st:st+127) = data(dataset(i):dataset(i)+127,1);
    st=st+128; 
end;
    
st=1;    
for i =1:length(dataset),
    plot(diff(time_seg(st:st+127))*24);
    pause;
    st=st+128; 
end;

    
 champ=load('C:\Manoj\projects\cej\stefan\EEJ_parameters_V3.0\EEJ_parameters_V3.0\CHAMP_scalar.txt');
inversion(data(:,1),data(:,5),10/24);
inversion(champ(:,1),champ(:,5),10/24);
c=inversion(champ(:,1),champ(:,5),10/24)
hh=logspace(1,10)
hh=logspace(0.1,1)
hh=logspace(0.1,1.1)
hh=logspace(0.1,2)
hh=logspace(0.5,2)
for i = 1:length(hh),
c(i,:)=inversion(champ(:,1),champ(:,5),hh(i)/24);
end;

%%

load  /data/backup/mnair/longp/data_limits 
champ=load('C:\Manoj\projects\cej\stefan\EEJ_parameters_V3.0\EEJ_parameters_V3.0\CHAMP_scalar.txt');
load  /data/backup/mnair/longp/eef_data_mod.mat eef;


eef1=eef;
eef1(:,6) = eef(:,6)-eef(:,15)/1e3; % remove climate 1e3 because the climate outputs in mV/m and eef(:,6) is V/m

data_index = eef_index;
N_seg=1;
N_data=1;
len = 128;
min_kp = 2;
clear data_sel time_sel
for i = 1: length(data_index),
%   champ_sel = champ(data_index(i,1):data_index(i,2),[1,6]);
    champ_sel = eef1(data_index(i,1):data_index(i,2),[1,6,9]);
    flag = 1;
    st=1;
    en=len;
    nn=1;
    
    while en <= length(champ_sel),
        timeax = diff(champ_sel(st:en,1))>0.07;

          
        if sum(timeax) > 0,
            fprintf('No %d %d %d %d\n',i,st,en,length(st:en))
            st= st+max(find(timeax==1))+1;
            en = st+(len-1);
            
        else,
         
                
            kp_sel = champ_sel(st:en,3);
            if nanmean(kp_sel) > min_kp,
            fprintf('Av Kp %5.2f %d %d %d\n',nanmean(kp_sel),st,en,length(st:en))
            data_sel(N_seg:N_seg+(len-1)) = champ_sel(st:en,2);
            time_sel(N_seg:N_seg+(len-1)) = champ_sel(st:en,1);
            time_marks(N_data,:) = champ_sel([st,en],1);
        
            
            N_seg = N_seg+len;
            N_data = N_data+1;

            end;

            st=st+len;
            en=en+len;
            nn=nn+1;
            
        end;
       
    end;
%     fprintf('i= %d nn = %d length = %d\n',i,nn,length(champ_sel));
end;

[pxx,f] = pwelch(data_sel*1e3,hanning(len),0,len,1/(92*60));
% loglog(1./(3600*f),pxx,'r');

%%
    %Notes enable some absolute indexing !!
%problem st en        
        

%resampling ACE data ACE data is 5 minutes intreval 
% resampling to median diff champ == 92.1974 minutes
%I chose 1/18 to have a defined X axis (1:18:end)
load  /data/backup/mnair/longp/OMNI_ELEC_new ace_all; % JUST ACE IEF

L = isnan(ace_all(:,2));
ace_all(L,2) = 0;
ace_down = resample(ace_all(:,2),1,18);
ace_inter = interp1(ace_all(1:18:end,1),ace_down,time_sel - 17/(60*24)); %17 min phase delay

[cxx,f] =mscohere(data_sel,ace_inter,hanning(len),0,len,1/(92*60));
semilogx(1./(3600*f),cxx,'r.-');
% 
% [cxx,f] =mscohere(data_sel(1:6400),ace_inter(1:6400),hanning(len),0,len,1/(92*60));
% semilogx(1./f,cxx,'r');
% 


% for i = 1:length(data_index),
%     
%     [trash1 ind1] = min(abs(eef(:,1) - champ(data_index(i,1),1)));
%         [trash2 ind2] = min(abs(eef(:,1) - champ(data_index(i,2),1)));
% %        fprintf('%d %d %d %d %5.4f %5.4f\n', data_index(i,1),ind1, data_index(i,2),ind2,trash1,trash2);
% eef_index(i,1) = ind1; eef_index(i,2)=ind2;
% end;
        


%%

load C:\Manoj\projects\ace\20080303\alldays JULI_SEG ACE_SEG;
JULI_SEG = JULI_SEG.*24.366*1e-3; %mV/m
[Txy1,F] = tfestimate(ACE_SEG,JULI_SEG,hanning(72),0,72,1/(5*60)); %Txy is agains calculated just to get
%the normalized F as required by the function invfreqz
%The TEST ACE data are corrected for 17 minutes
%[b,a]=invfreqz(Txy,F,37,0);
[b,a]=invfreqz(Txy1,F,7,4,[],3000); %better

load c:\manoj\projects\ace\OMNI_ELEC ace_all; % JUST ACE IEF
load c:\manoj\projects\longp\selected_eef_champ time_sel data_sel len time_marks N_data; %time_mark contains the start and end of each segment
N_seg=1;
time_mark = time_marks-17/(60*24); %17 minute phase delay
for i = 1:N_data-1,
    [trash1 ind1] = min(abs(ace_all(:,1) - (time_mark(i,1)- 6/24))); %6/24 to get some time b4 
    [trash2 ind2] = min(abs(ace_all(:,1) - time_mark(i,2)));  
    ace = ace_all(ind1:ind2,2);
    L = isnan(ace);
    ace(L) = 0;
    ace_time = ace_all(ind1:ind2,1);
    pred_julia = filter(b,a,ace)*2; %assumption ace raw data are 5 minutes interval ?!!%$
    fprintf('%d %d\n',i,length(ace));
    plot(ace_time,pred_julia);
    hold on;
    plot(time_sel(N_seg:N_seg+(len-1)),(data_sel(N_seg:N_seg+(len-1))-nanmean(data_sel(N_seg:N_seg+(len-1))))*1e3,'r');%1e3 to make it mV/m
    pred_ace_at_champ = interp1(ace_time,pred_julia,time_sel(N_seg:N_seg+(len-1)));
    

    
xx=corrcoef(pred_ace_at_champ(1:end-1),(data_sel(N_seg:N_seg+(len-2))-nanmean(data_sel(N_seg:N_seg+(len-2))))*1e3);
title(sprintf('%5.3f',xx(2,1)));
    
    pause;hold off;    
    N_seg = N_seg+len;

end;
    
%%
uiopen('C:\Manoj\projects\longp\julia_champ_spectra_new.fig',1);
uiopen('C:\Manoj\projects\longp\JULIA_CHAMP_ACE_COH.fig',1);
uiopen('C:\Manoj\projects\longp\Julia_Ace_Champ_tf.fig',1);    
    
%the following script (upto the next comment) is to interpolate CHAMP data to equal interval

if ~exist('eef','var'),
    display('loading EEF data');
load c:\manoj\projects\longp\eef_data_mod eef;
eef1=eef;
eef1(:,6) = eef(:,6)- eef(:,15)/1e3; % remove climate 1e3 because the climate outputs in mV/m and eef(:,6) is V/m
end;

des_int = 2;%desired sampling interval in hours 
len = 128;%window length for fft, coherence or transfer function estimate
cl_c = 'r.-';
win_dow = parzenwin(len);
overlap = 1;

%bartlett, blackman, chebwin, hamming, hann, kaiser,parzenwin

%des_int ==2 and len == 16 gave best result !!!!! parzen window or bartlett
%gave good res coh and tf but these windiws ive somewhat jerkey e-field
%(JULIA-CHAMP). No overlap resulted in coh and tf values somewhat reducing
%- no change to eef spectra

%time_axis_desired = eef1(1,1):1/24:eef1(end,1);
  
time_axis_desired = eef1(1,1):des_int/24:eef1(end,1);
eef_at_equal_interval = interp1(eef1(:,1),eef1(:,6),time_axis_desired);
%eef_down_to_2 = resample(eef_at_equal_interval,1,2);

[pxx,f] = pwelch(eef_at_equal_interval*1e3,win_dow,overlap,len,1/(des_int*3600));
figure(1);
loglog(1./(3600*f(2:end-1)),pxx(2:end-1),cl_c,'LineWidth',2);
hold on;


if ~exist('ace_all','var'),
    display('loading ACE data');
load c:\manoj\projects\ace\OMNI_ELEC ace_all; % JUST ACE IEF

L = isnan(ace_all(:,2));
ace_all(L,2) = 0;
end;

%ace_all(:,2) = normrnd(nanmean(ace_all(:,2)),std(ace_all(:,2)),[1,length(ace_all)]);
% The above line was to test the "null hypothesis". A random data against
% the CHAMP. Found that there is no coherence between them (<0.002)

ace_down = resample(ace_all(:,2),1,des_int*60/5);
%ace_inter = interp1(ace_all(1:18:end,1),ace_down,time_axis_desired - 0/(60*24)); %17 min phase delay
ace_inter = interp1(ace_all(1:des_int*60/5:end,1),ace_down,time_axis_desired - 17/(60*24)); %17 min phase delay

[cxx,f] =mscohere(eef_at_equal_interval,ace_inter,win_dow,overlap,len,1/(des_int*3600));
figure(2);hold on;
semilogx(1./(3600*f),cxx,cl_c);




%JULI_SEG = JULI_SEG.*24.366*1e-3; %mV/m
[Txy,F] = tfestimate(ace_inter,eef_at_equal_interval*1e3,win_dow,overlap,len); %Txy is agains calculated just to get
%the normalized F as required by the function invfreqz
%The TEST ACE data are corrected for 17 minutes
%[b,a]=invfreqz(Txy,F,37,0);
%[b,a]=invfreqz(Txy1,F,7,4,[],3000); %better

figure(3);hold on;
loglog(1./(3600*f),abs(Txy),cl_c)


%%

%load c:\manoj\projects\longp\CHAMP_ACE_DATA_F1 eef_at_equal_interval ace_inter len des_int win_dow overlap;
load  /data/backup/mnair/longp/CHAMP_ACE_DATA_F1 eef_at_equal_interval ace_inter len des_int win_dow overlap;
[Txy_long,F_long] = tfestimate(ace_inter,eef_at_equal_interval*1e3,win_dow,overlap,len,1/(des_int*3600)); %Txy is agains calculated just to get


load  /data/backup/mnair/longp/alldays JULI_SEG ACE_SEG;
JULI_SEG = JULI_SEG.*24.366*1e-3; %mV/m
[Txy_short,F_short] = tfestimate(ACE_SEG,JULI_SEG,hanning(72),0,72,1/(5*60)); %Txy is agains calculated just to get
%the normalized F as required by the function invfreqz
%The TEST ACE data are corrected for 17 minutes
%[b,a]=invfreqz(Txy,F,37,0);

Txy_short(1:2) = [];
F_short(1:2) = [];
%this eleminates the logest period (6 hours) and 0 frequencies from
%JULI-ACE Tf



[F_all,IA,IB] = union(F_long,F_short);
[A,B] = sort([F_long(IA);F_short(IB)]);

%A and F_all should be the same;

T_all = [Txy_long(IA);Txy_short(IB)];
T_all = T_all(B); %sort the Txy according the the frequency mix;
F_scaled = pi*(F_all./max(F_all));

T_all(1) = 0; % removes the DC average from the transfer function...

%
%Wt = (abs(T_all))./max(abs(T_all)); %another way would be to use Coh

Wt = zeros(size(T_all));
Wt([2,11,42]) = 1;% The values corresponds to start, middle (2 hours) and end of frequencies
[b_n,a_n]=invfreqz(T_all,F_scaled,4,4,Wt,3000,0.001); %4 and 4 found out from trial and error



%%
% The following script was used to find optimum numbers of coeffieicits ( b
% and a) for filter. Also tested the weighting algorithm etc.
%I tried to fit the frequeiciens around 2 hours 
%as well as the highest and lowest frequencies while desiging the filter
%for this the weight Wt vector was used. The values at 2,11 ans 42 were set
%to 1 and the rest to zeros. It may be possible to have values ranging from
%0-1, but it is beyond the scope of this work.
%In this set-up, an=4.bn=4 resulted in best estimate

len_win = 360;
step = zeros([1,len_win]);
step(13:end) = 1; % positive step

for i = 1:5,
    for j = 1:5,
        L = Wt > 0;
        [b,a]=invfreqz(Txy_short,pi*(F_short./max(F_short)),7,4,[],3000); %better
        [b_n,a_n]=invfreqz(T_all,F_scaled,i,j,Wt,3000,0.001); %better
        subplot(211);
        plot((1:5:len_win*5)./60,filter(b,a,step)*1.5,'b.-','LineWidth',1);
        hold on;
        plot((1:5:len_win*5)./60,filter(b_n,a_n,step),'k.-','LineWidth',1);
        hold off;
        [h,w] = freqz(b,a,F_all,1/(5*60));
        [h_n] = freqz(b_n,a_n,F_all,1/(5*60));
        
        subplot(212);
%         semilogx(F_all,real(h),'b',F_all,imag(h),'b');hold on; 
%         semilogx(F_all,real(h_n),'k',F_all,imag(h_n),'k')    ;
%         semilogx(F_all,real(T_all),'r.',F_all,imag(T_all),'ro');
%         semilogx(F_all(L),real(T_all(L)),'r*',F_all(L),imag(T_all(L)),'r*');

        loglog(F_all,abs(h),'b');hold on; 
        loglog(F_all,abs(h_n),'k');   ;
        loglog(F_all,abs(T_all),'r.');
        loglog(F_all(L),abs(T_all(L)),'r*');

%         semilogx(F_all,real(h),'b');hold on; 
%         semilogx(F_all,real(h_n),'k');   ;
%         semilogx(F_all,real(T_all),'r.');
%         semilogx(F_all(L),real(T_all(L)),'r*');

        
        hold off;
        title(sprintf('an = %d bn = %d',i,j));
        fprintf('%d %d\n',i,j);
        pause;
    end;
end;

% 
% 
% 
% 
% for i = 1:20,
%     for j = 1:20,
%         
%         [b,a]=invfreqz(T_all,F_scaled,i,j,[],30); %better
%         [b_n,a_n]=invfreqz(T_all,F_scaled,i,j,Wt,30); %better
%         [h,w] = freqz(b,a,F_all,1/(5*60));
%         [h_n,w_n] = freqz(b_n,a_n,F_all,1/(5*60));
%         err1 = sqrt(mean(abs(h.*conj(h)-T_all.*conj(T_all))));
%         err2 = sqrt(mean(abs(h(7:17).*conj(h(7:17))-T_all(7:17).*conj(T_all(7:17)))));
%         fprintf('i= %d j = %d ERR = %f  %f\n',i,j,err1,err2);
%         error_mt(i,j) = err1;
%         error_mt1(i,j) = err2;
%     end;
% end;
% 
% 
% 



%%




%interpolate CHAMP EEF to a constant sampling interval, respecting the data
%gaps. A way to do this would be to make a time axis with desired interval,
%delete all the points which are away from a thrashold of the observed
%data's time axis. For example, let 1789.1 1789.2 ... 1789.9 is the desired time axis
%if as observation was made at 1789.48, 1789.62 the desired time scale's
%all other entries other than 1789.5 and 1789.6 will be deleted.. Then a
%bilinear interplotaion of observed data to these points (1798.5 and
%1789.6) will be performed. Thrashold will be time interval (or half ?) of
%the desired time axis


load c:\manoj\projects\longp\eef_data_mod eef;
eef1=eef;
eef1(:,6) = eef(:,6)-eef(:,15)/1e3; % remove climate 1e3 because the climate outputs in mV/m and eef(:,6) is V/m
time_axis_desired = eef1(1,1):1/24:eef1(end,1);
eef_at_equal_interval = interp1(eef1(:,1),eef1(:,6),time_axis_desired);



for i = 1:length(time_axis_desired),
    
[trash2 ind2] = min(abs(eef(:,1) - time_axis_desired(i)));  
deviation_time_axis(i,1) = trash2;
deviation_time_axis(i,2) = ind2;

end;
    

L = deviation_time_axis(:,1) > 1/24;
time_axis_desired(L) = [];
eef_at_equal_interval(L) = [];

eef_at_equal_interval = interp1(eef1(:,1),eef1(:,6),time_axis_desired);

load c:\manoj\projects\longp\data_limits;
data_index = eef_index;



for i = 1:length(data_index),
    
    [trash1 ind1] = min(abs(time_axis_desired - eef1(data_index(i,1),1)));  
    [trash2 ind2] = min(abs(time_axis_desired - eef1(data_index(i,2),1)));  
    
    new_data_index(i,:) = [ind1,ind2];
end;


%after the time axis is equal spaced and the no-data-poits are deleted, let
%us make the time and data into pieces of equal number of continous time
%series segments

N_seg=1;
N_data=1;
len = 32;
min_kp = 2;
clear data_sel time_sel

%%
for i = 1: length(data_index),
    clear champ_sel;
    champ_sel(:,2) = eef_at_equal_interval(new_data_index(i,1):new_data_index(i,2));
    champ_sel(:,1) =     time_axis_desired(new_data_index(i,1):new_data_index(i,2));
    flag = 1;
    st=1;
    en=len;
    nn=1;
    
    while en <= length(champ_sel),
        display('here');
        timeax = diff(champ_sel(st:en,1)) > 0.0417; %(1/24)

          
        if sum(timeax) > 0,
            fprintf('No %d %d %d %d\n',i,st,en,length(st:en))
            st= st+max(find(timeax==1))+1;
            en = st+(len-1);
            display('herenthere');
        else,
         
                
%             kp_sel = champ_sel(st:en,3);
%             if nanmean(kp_sel) > min_kp,
            fprintf(' %d %d %d\n',st,en,length(st:en));
                    display('there');
            data_sel(N_seg:N_seg+(len-1)) = champ_sel(st:en,2);
            time_sel(N_seg:N_seg+(len-1)) = champ_sel(st:en,1);
          %  time_marks(N_data,:) = champ_sel([st,en],1);
        
            
            N_seg = N_seg+len;
            N_data = N_data+1;

    
            st=st+len;
            en=en+len;
            nn=nn+1;
            
        end;
       
    end;
%     fprintf('i= %d nn = %d length = %d\n',i,nn,length(champ_sel));
end;
[pxx,f] = pwelch(data_sel*1e3,hanning(len),0,len,1/(60*60));
loglog(1./(3600*f),pxx,'r');



% The above processing resulted low-biased estimate of  high frequency part of the spectra 
%%


%% March 12, 2012 
% estimating the error of the transfer functions
% Following the paper by MacMynowski and Tziperman (Article submitted to
% Royal Society)


% find the coherence between the signals
% uiopen(' /data/backup/mnair/longp/tf_db_magnitue_julia_champ.fig',1);
hold on;
load  /data/backup/mnair/longp/CHAMP_ACE_DATA_F1 eef_at_equal_interval ace_inter len des_int win_dow overlap;
[Txy_long,F_long] = tfestimate(ace_inter,eef_at_equal_interval*1e3,win_dow,overlap,len,1/(des_int*3600)); %Txy is agains calculated just to get
[cxx,f] =mscohere(ace_inter,eef_at_equal_interval*1e3,win_dow,overlap,len,1/(des_int*3600));
Err = sqrt( 1/(2*(length(eef_at_equal_interval)/len)) .* ( (1-cxx)./cxx ) ) .* abs(Txy_long);
Err_long_lg = 0.434 * Err ./ abs (Txy_long);
errorbar(log10((1./F_long(2:end))./(3600)), log10(abs(Txy_long(2:end))),Err_long_lg(2:end)/2,Err_long_lg(2:end)/2,'r');
hold on;
%%
load  /data/backup/mnair/longp/alldays JULI_SEG ACE_SEG;
JULI_SEG = JULI_SEG.*24.366*1e-3; %mV/m
[Txy_short,F_short] = tfestimate(ACE_SEG,JULI_SEG,hanning(72),0,72,1/(5*60)); %Txy is agains calculated just to get
[cxy_short,f_short] = mscohere(ACE_SEG,JULI_SEG,hanning(72),0,72,1/(5*60));

Err_short = sqrt( 1/(2 * 265) .* ( (1-cxy_short)./cxy_short ) ) .* abs(Txy_short);
Err_short_lg = 0.434 * Err_short ./ abs (Txy_short);

errorbar(log10((1./F_short(2:end))./(3600)), log10(abs(Txy_short(2:end))),Err_short_lg(2:end),Err_short_lg(2:end),'b');

%% March 17, 2012
% Fix issues with the CHAMP eef data gaps


% the climatology is limted to > 7 and < 17 LT
load  /data/backup/mnair/longp/eef_data_mod.mat eef;
tol = 1.8/24; % max time interval tolerated
min_time_length = 8/24; % minimum time length required in decimal days
np = 0;
lt_start = 7;
lt_end = 17;
nd = 1;

for i = 1: 36665,
    
    if (eef(i+1,1) - eef(i,1) <= tol && eef(i,4) > lt_start && eef(i,4) < lt_end ...
            && eef(i,1) - eef(i - np ,1) <= min_time_length )
        np = np + 1;
    else
       % if (np >=16)
        if ( eef(i,1) - eef(i - np ,1) >= min_time_length )
            data_index(nd,2) = i;
            data_index(nd,1) = i - np;
            np = 0;
            nd = nd + 1;
        else
            np=0;
            
        end;
    end;
    
end;
save  /data/backup/mnair/longp/eef_data_mod.mat eef data_index min_time_length;

%% Further dviding the data index to smaller pieces
% NOT used anymore (the above script does this work too)
% required_length = 32;
% required_time_length = 32/24; % required time length
% flag = 1;
% nd = 1;
% 
% for i = 1 : length(data_index)
%     
%     nseg = floor( ( eef(data_index(i,2),1) - eef(data_index(i,1),1) ) / required_time_length ) ;
%     
%     for j = 1 : nseg,
%         data_index_new(nd,1) = data_index ( i , 1 ) + ( j - 1) * required_length  ;
%         data_index_new(nd,2) = data_index ( i , 1 ) +  j * required_length - 1;
%         nd = nd + 1;
%     end;
% end;
%         
%     
    
   

%% Filtering Jicamarca ISR data in the same way as the CHAMP EEF
% data. Idea is to organize the data into continous samples of a given
% lnegth 
load ([data_dir 'ace/jcamarca_isr_day_night.mat']);
tol = (5/60)/24 + 0.0001; % max time interval tolerated
min_time_length = 32/24; % minimum time length required in decimal days
np = 0;

nd = 1;

for i = 1: 25464,
    
    if (fday(i+1) - fday(i) <= tol & fday(i) - fday(i - np ) <= min_time_length )
        np = np + 1;
    else
       % if (np >=16)
        if ( fday(i) - fday(i - np) >= min_time_length )
            data_index(nd,2) = i;
            data_index(nd,1) = i - np;
            np = 0;
            nd = nd + 1;
        else
            np=0;
            
        end;
    end;
    
end;
save  /data/backup/mnair/longp/jicamarca_isr_data.mat fday drift data_index min_time_length;